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Abstract 

We analyze the complexity of the quantum optimization algorithm based on adiabatic evolution 
for the NP-complete set partition problem. We introduce a cost function defined on a logarithmic 
çf} [ scale of the partition residues so that the total number of vàlues of the cost function is of the order 

of the problem size. We simulate the behavior of the algorithm by numerical solution of the time- 
dependent Schròdinger equation as well as the stationary equation for the adiabatic eigenvalues. 
The numerical results for the time-dependent quantum evolution indicate that the complexity of the 
algorithm seales exponentially with the problem size. This result appears to contradict the recent 
numerical results for complexity of quantum adiabatic algorithm applied to a different NP-complete 
problem (Farhi et al, Science 292, p.472 (2001)). 

Oh! 

c3 ; 1 Introduction 

cr. 

Most common computationally intensive tasks encountered in practice may be formulated as combi- 
^ ■ natorial optimization problems (COPs), many of which are found to belong to the algorithmic class 
nondeterministic-polynomial complete (NP-complete) |ÏJ. The NP-complete problems are computa- 
tionally hard - they are characterized (in the worst cases) by exponential sealing of the running time or 
memory requirements with the problem size. A special property of the class is that any NP-complete 
problem can be converted into any other NP-complete problem in polynomial time on a classical com- 
puter; therefore, it is sufficient to find a deterministic algorithm that can be guaranteed to solve all 
instances of just one of the NP-complete problems within a polynomial time bound. 

An instance of a COP of size n may be encoded using bit strings z = zq Z\ ■ ■ ■ z n —it z j = 0, 1, 
with a corresponding value of the cost function (or "energy") E = E z for each string. The objective 
is to find the bit string(s) with the minimum cost (and the corresponding cost value). In quantum 
computation, bits Zj are replaced by spin-^ qubits; the qubit states |0) and |1) are eigenstates of the 
■izz component of the j'-th spin, respectively. The Hilbert space of a quantum register with n qubits is 
spanned by N = 2 n basis vectors \z) = \zq) (g) • • • (g) \z n -\). 



2 Optimization by Adiabatic Quantum Evolution 

Following ||, U, we consider a quantum evolution of duration T based on the time-dependent 
Hamiltonian 

H(t) = a(s)V + (3(s)H P , s = s{t). (1) 

Here, Hp = ^2 z E z \z)(z\ is the "problem" Hamiltonian that embodies the problem structure in its 
energy spectrum and eigenstates, the summation being performed over all ./V ro-bit strings, and V is a 
"driver" Hamiltonian that is constructed in such a way as to cause transitions between those states - 
essentially an Ising-type spin Hamiltonian corresponding to 1- and 2-gate operations: 

re— 1 -, re— 1 



^Bioi-l È -VX- (2) 



i=0 i,j=0 

Coeficients a (s (i)) and f3(s(t)) vary in time in such a way that at the initial instant of time 
H(0) = V and at the final instant H(T) = Hp. A particular choice of the coefficients is ||] 

a(s) = l-s, /3(s) = s, s(t)=t/T (3) 

The total Hamiltonian (|^) produces a nontrivial quantum evolution from some initial (superposition) 
state V'(O) to a final (solution) state ip(T). If no knowledge about the solution is available a priori, 
then the initial state may be chosen as the symmetric state (cf. [||. ) 

2 n -l 

V-(0) = 2~«/ 2 ^ |z>. (4) 

z=0 

This choice is appropriate provided (Q) is a ground state of V (e.g., Bi, Jy > 0). Now, if T is sufficiently 
large, then functions a(t/T) and f3{t/T) vary in time slowly and the system will remain in the instan- 
taneous (adiabatic) ground state of H(t) during its entire evolution < t < T (cf. Q). Accordingly, 
ip(T) will be a superposition of states \z) corresponding to the ground state of Hp. It is clear that in 
this case a measurement performed on the quantum register at t = T will find with certainty one of the 
solutions of COP. In this case the complexity of the quantum algorithm is determined by its duration 
T. If we expand the wavefunction of the system ip(t) in the basis of the adiabatic eigenfunctions ^(í) 
of the Hamiltonian H(t) 

H(t)V k (t)=g k (t)ï> k (t), fe = 0,1,... ,2 n -l. (5) 

2 n -i / ,í N 

(6) 



m= Y\ C k (t)V k (t)exp(-i/n f dt'g k (t'\ 



(7) 



then adiabatic approximation corresponds to ip(t) oc ^o(t) ( U P 1° the oscillating phase factor). Coeffi- 
cients Cfc(í) with k > correspond to nonadiabatic corrections. Using perturbation theory in the basis 
of eigenfunctions í'fc(í) the total probability p n - a d{t) of not finding the system at the instant t in its 
adiabatic ground state equals 

Here we used the explicit form of coeficients a, (3 given in (|3|). It is seen that during the quantum 
evolution coefncients C&(í) ~ s/s ~ 1/T and the largest admixture of the exited states into the total 
superposition occurs at the instant of time when one of the exited levels closely approaches the 
ground state (avoided-crossing). From here the overall criterion for the adiabatic evolution can be 
expressed in the well-known form 

V = ^~«h (9) 

where Ag m i n is the closest approach of the ground state to one of the excited states during the evolution 
- a minimum gap- and V is the characteristic energy scale for the matrix elements of V. We note that 
although instantaneous nonadiabatic corrections (||) are quadratic in the parameter rj (near the avoided 
crossing) the probability of nonadiabatic transitions W n - a d away from the ground state is exponentially 
small in rj B. This probability is defined on an infmite time axis and its logarithm is proportional to 
the imaginary part of the integral along the contour in the complex plane of t that begins and ends on 
the real time axis and loops around the complex branching point t* 



W n . ad exp I -|Im^> g (t)dt\j . (10) 
Here t* correspond to one of the roots of the equation 

9k(t*) = so (O (11) 



that provides the smallest value for the exponential in ( |Ï0|) (out of all possible complex solutions of 
( p^ ) for different excited states ^?k)- I n a Standard (Landau-Zener) theory of nonadiabatic transitions 
the value of the exponent is approximately of the order of the parameter 77 in (|9[) , and therefore it is the 
size of the minimum gap Ag m i n that determines the condition for X 1 and hence the complexity of the 
quantum adiabatic search algorithm according to [§]. We note finally that, as pointed out in §, the 
improved complexity of the adiabatic algorithm is determined by the instantaneous rate s(t)/s(t) of 
the variation of the control parameter s(t) near the avoided crossing. We will not discuss in this paper 
such modifications and focus primarily on intrinsic properties of the quantum system in question. 
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3 Set Partition Problem 



In this paper, we will analyze the complexity of the adiabatic quantum optimization for the set partition 
problem (SPP), which is one of the bàsic NP-complete problems of theoretical computer science Q. 
The optimization version of SPP is to partition a set of n positive integers {ao,a\, . . . ,a n _i} into 
two disjoint subsets A\ and Ai such that the "residue" | Y^a-eAi a i ~ ^L·ia^Ai a i\ * s m i n i m ized. The 
complexity of the problem substantially depends on the size of the integers ctj (see below). It is often 
customary for the analysis of the random instances of the problem to introduce finite-precision rational 
numbers Oj that are independently and identically distributed (i. i. d.) in the unit interval (0, 1]. 

ctj < 2 b Vj, dj = 2~ b aj £ (0, 1] (12) 

Here b is the total number of bits used to represent the numbers a,j. The vàlues of the partition can 
be encoded in binàries by attaching "sign" bits Sj to the numbers a,j. The partition residue can be 
defined as |Í2 Z | where 

n-l 

íl z = s j a ji s j = 1 ~~ 2zj = ±1 (zj = 0, 1) (13) 

3=0 

Here VL Z is a signed partition residue. We note that by definition the problem is symmetric: two bit 
strings that can be obtained from each other by flipping all the bits (sj — > —Sj) correspond to two 
vàlues of Í2 Z that differ only in sign. We note that the minimum-residue partition(s) may be thought 
of as the ground state(s) of the following spin Hamiltonian |7| 

n-l 

íl z = aidjSiSj. (14) 
i, 3=0 

This is an infinite range Ising spin glass with Mattis type antiferromagnetic coupling, Jj,- = —a^aj. 
Infinite range coupling clearly represents a major problem with direct ('analog') physical implemen- 
tation of this Hamiltonian on a quantum computer. Therefore one can consider using an oracle-type 
cost function E(z) = |íí z | to implement the problem Hamiltonian in ([]]) for SPP. The corresponding 
unitary transformation will multiply the basis states \z) by phase factors exp (— iAt E(z)) during the 
elementary discrete steps of the 'continuous-time' adabatic quantum optimization (^). Although this 
approach is natural for the satisfiability problem Q it has a serious limitation for SPP (as well as some 
other NP-complete problems like integer programming, where the precision of integers is of central 
importance). To demonstrate this point we need to consider the density of states of the partition 
residues. 
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3.0.1 Density of states 

We define the density of states for a given instance of SPP as follows 

p (o) = ^2ó(n-n z ). (15) 

z 

The exact form of p(£l) depends on a given instance of SPP (i.e., a particular set of numbers dj). 
However we introduce a coarse-grained density of states 



p(0) = ^/ ,>kkk (1.0) 



cCl+AÜ 
IQ 

where averaging is over an interval of AQ whose size will be determined below. Using (^) and (15) 
we can rewrite this expression in the form 

m = — / : I(w) dw, I(w) = H cos(a jW ). (17) 

Z7T J -OO lW j_q 

Note that I(irk2 b ) = (— l) k , k = 0, ±1, . . . and I(w) has very sharp màxima (minima) at those points. 
In their vicinities the integral in ( |Ï7l) can be evaluated by steepest descent method for any given problem 
instance. The sum over the contributions from different saddle points was obtained by Mertens in 
his derivation of the partition function for the corresponding spin glass model. We emphasize however 
that I(w) can have múltiple sharp resonances at the intermediate points \w\ < 2~ b . The positions of 
these resonances are at the múltiples of n/q where q is an approximate greatest common divisor (g.c.d.) 
of the set of n numbers a,j such that a,j = fjq + rj where fj are integers and rj are residues of the 
division. Provided that most of the residues are sufhciently small 

2 



7T v ^ 2 



the function I(w) will have steep peaks at those points. It can be shown that in the general case the 
value of the approximate g.c.d. for a set of n b — bit numbers inside the unit interval scales as 2~ n for 
n < b. Obviously it equals 2~ b for n > b. In what follows we will be interested in the high-precision 
case n < b. We choose the size of the averaging window Aíi ^> 2~ n and this introduces a cut-off in the 
integral (|Ï4|) at 

\w\ <C tt/AÍÍ <C 2 n 

It follows from above that in this case the vàlues of the g.c.d. will lie outside the cutoff and corresponding 
resonances will not contribute to the integral. The value of the integral can be estimated near the single 
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remaining maximum at w = 0. The width of the maximum near that point is ów n l l 2 <C 1 and 
therefore the window function in (|ï~3|) works as a step function in that region. Finally we obtain 



2 n 
y/2irna 2 



/ 2 \ i «- 1 



i=o 

Here the variance cr is a "self-averaging" quantity, and the coarse-graining is performed over an interval 
much larger than the characteristic separation between neighboring partition residue vàlues 

AE ~ y/n 2~ n (19) 

yet much smaller than the scale of variation of p(í7) : AE <C Aíl <C y/n. We note that in the high- 
precision regime (n < b), partition residues íí z are irregularly spaced and well separated from each other 
(on the scale of 2~ b ). However this structure is being averaged out in (|ï~8| ) and the result indicates 
that, in general, no more structure exists on a scale 3> 2~ n other than that given by the Gaussian 
distribution in (18). We note that this distribution is usually obtained for the SPP using averaging 



over different instances of the problem (cf. j|, @(&)); here we recovered it as a coarse-grained 
distribution for a given instance which is more consistent with our goal of studying the complexity of 
the adiabatic quantum optimization algorithm |Ï0| ]. 

3.1 Cost function 

Computational complexity of SPP depends critically on the number of bits ò: numerical simulations 



with independent identically distributed (i. i. d.) random ò-bit numbers aj show [11, [Ï2|] that the solution 
time grows exponentially with n for 2 n ~ b <C 1 (high-precision, computationally 'hard phase'), and 
polynomially for 2 n ~ b 3> 1 (low-precision, computationally 'easy phase'), exhibiting a behavior similar 
to a phase transition 0|(a). 

In the low-precision phase, vàlues of Q z are equally spaced (in 2~ b ) and strongly degenerate each 
corresponding to (roughly) 2 b ~ n 3> 1 number of bit-strings. This degeneracy grows exponentially with 
n if b remains fixed. The total number of solutions with zero residues accumulate correspondingly and 
this is why the complexity eventually becomes polynomial in n. The quantum algorithm suggested in 



[13] directly computes the density of states ( |Ï5| ) of the SPP and is emcient in finding the number of 
solutions in the low-precision case. In this case it is also feasible to use a cost function E(z) = |fi z | 
(provided the number of possible vàlues does not grow exponentially with n). 

The situation is qualitatively different in the high-precision case. Implementation of the approach 
based on the above cost function will require a quantum computer using exponentially high precision 
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physical parameters (external fields, etc) to control small differences in the phases of unitary trans- 
formations on the scale at least 2~ n . This is a technical difference from the constraint satisfaction 
problem in which the cost function generally takes only the set of vàlues that scales polynomially with 
n; the size of the set equals the total number of constraints m (the computationally most difficult 
case corresponds to m ~ n and the case of m ~ 2 n is not of general interest there). To avoid the 
above restriction in the implementation of the adabatic quantum optimization algorithm for SPP, we 
introduce a cost function E(Q) based on a logarithmic scale of the partition residue vàlues: 

E(Q) = 0, for < |0|/A < 1, (20) 
E{Ü) = k, for 2 fc ~ 1 < |fi|/A < 2 k , k = 0, 1, . . . L, 

n-l 

2 L ~ 1 < A/A < 2 L , A = ^2dj. 

j=0 

Since the density of states is linear at fi <C \fn number of states dk per energy level E(Q Z ) = k will 
grow exponentially with k in that range (~ 2 k ). The total number of levels L depends on the value of 
A. Using the density of states (|Ï£|) for Í2 Z <C \fn one can set 

A = Vn~2- n K. (21) 

where K is some fixed number (a few dozen) independent of n. The number of ground states of the 
problem Hamiltonian 

L-l 

H P = Y J E{ïl z )\z){z\ (22) 

k=0 

approximately equals K, and the total number of energy levels L is close to n. It can be estimated 
from ( |2"0| ) that n — L ~ log2K. The distribution of the low-lying states dk with cost function ( |20| ) is 
somewhat similar to that in the slightly underconstrained cases of the satisfiability problem. 



4 Results 

To study the complexity of the adiabatic quantum optimization algorithm for SPP, we numerically 
integrate the time dependent Schròdinger equation with the Hamiltonian H(t) (|^), (0) in which we 
set Bi = 1 and Jij = 0. We start from the symmetric initial state (||) and integrate the Schròdinger 
equation in the interval < t < T. Unlike the approach adopted in |J we do not set a-priori a value 
of success probability. Instead we introduce a complexity mètric for the algorithm 

Cm - ,23) 
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Figure 1: C(T) vs T for n=15, precision b=25 bits. Point 1 on the figure corresponds to the minimal 
value of complexity; the corresponding vàlues are T* = 22.67, po(T*) = 0.15 and do = 22. At Point 2 
the total population of the ground level has aheady reached po(T) = 70%. 

Here Po(T) is the total probability of finding the system in its ground level (with E(íl z ) = 0) at the 
end of the algorithm, t = T, and d® is the number of states at the ground level. The algorithm has to 
be repeated on average do/po(T) number of times to reach success probability ps 1. A typical plot of 
C(T) for an instance of SPP with n=15 numbers is shown in Fig. |ï[ At very small T the wavefunction 
is close to the symmetric initial state and the complexity is ~ 2 n . The extremely sharp decrease in 
C(T) with T is due to the buildup of the population po(T) in the ground level as quantum evolution 
approaches adiabatic limit. At certain t = T* the function C(T) goes through the minimum: for 
T > T* the decrease in the number of trials do/po(T) does not compensate anymore for the overall 
increase in the runtime T for each trial. The minimal complexity C* = C(T*) is defined via one 



dimensional minimization over T for a given problem instance [14|. In Fig.g we plotted the data for 
optimal complexities C* at different vàlues of n on a logarithmic scale. Vertical sets of points on the plot 
indicate the results for all simulation data we currently have for each n. The results indicate that the 
median value of complexity C* scales exponentially with n; linear fit to the graph gives logC* ~ 0.56 n. 
This corresponds to the scaling law C* ~ 2 ' 8 ™. The exponential behavior of the algorithm clearly 
manifests itself for the larger vàlues of n > 11. The scatter in the vàlues of logC* appears to decreases 
with n however this result is probably due to the smaller number of data points available for larger n 
vàlues. 

In Fig. p| we show the distribution of the probabilities |(z|?/'(r)| 2 for different vàlues of T for an 
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Figure 2: C* vs n, precision b=25 bits. Percentage figures correspond to the total population at the 
ground level at minimal complexity for a given n (e.g. 25% for n=9). Numbers below the vertical sets 
of points for each n show the number of trials (e.g., 97 for n=8, 6 for n=16, etc). Numbers on the top 
indicate the average vàlues of do for all trials at a given n. Median vàlues for the complexity for each n 
are shown with red squares. The line is the least squares fit of an exponential function to the median 
vàlues between n=ll and n=17. 



9 




Figure 3: | (z\tp(T)\ 2 vs z for one instance of SPP with n=15. Note that the vàlues of the index number 
on the horizontal axis correspond to the positions of different bit-strings z sorted with respect to the 
partition residue vàlues |íí 2 | (in increasing order). Index number corresponds to the smallest partition 
residue. The number of states at the ground level is cLq=22. Curve shown in red corresponds to the 
value of T = T* = 32 (minimal complexity). Curve shown in green corresponds to the value of T=8 
and black color curve corresponds to T=90. 

instance of SPP with n=15 (plots for different T shown with different colors)and precision b=25 bits. 
Vàlues of z are ordered with respect to the corresponding vàlues of the partition residues |í) z |. It is 
clearly seen on logarithmic scale that probability distribution forms 'steps' corresponding to different 
vàlues of the cost function E(íl z ) defined in (p0|). Within each step, the distribution of probabilities 
does not reveal any structure. The same property holds also for intermediate times (i < T). Detailed 



analytical results [15] indicate that it is this absence of structure in ip z (t) that is responsible for the 



exponential complexity of the algorithm. 
The Stationary Schròdinger equation 

In addition to solving the time-dependent Schròdinger equation we also analyzed the adiabatic solutions 
of the stationary Schròdinger equation with the same form of the Hamiltonian H(t) (|ï|) as above. Our 
preliminary results were obtained using Mathematica for modest vàlues of n < 10. The results for 
n=10 are shown in Figs. |] and ||. Figure |B| represents the magnified part of Fig. |] near the avoided 
crossing region. Adiabatic eigenvalues were computed for different vàlues of the scaled time parameter 
s = t/T G (0, 1). The sòlid line represents the evolution of the ground state eigenvalue between s = 
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and s = 1. The vertical sets of points correspond to excited adiabatic levels for a given s. At the 
beginning (s = 0) eigenvalues correspond to those of the Hamiltonian V: equally spaced levels -n,-n+2, 
. . . , n, corresponding to different number of spin excitations along the x quantization axis. The first 
excited state is n— fold degenerate, the second is n{n — 1), the k-th exited state is (^)-fold degenerate, 
etc. For s > the degeneracy is removed. For s — > 1 the eigenvalue spectrum is the one for the 
problem Hamiltonian Hp (in [22]) . We have shifted the energy reference in the Hamiltonian Hp by — n 
(cf. also ([ï|)) to match the energy scale for the symmetric case which emphasizes the avoided crossing 
region. In our case the ground state was 13-fold degenerate and the corresponding eigenvalues merge 
at s — > 1. The close approach of these eigenvalues is not relevant for the minimum-gap analysis since 
they all end up in the same final level. However the minimum separation of the instantaneous adiabatic 
ground state eigenvalue from the excited state eigenvalues that do not end up on the same ground level 
at s = 1 is clearly seen in the figures. Note that the size of this separation is much greater than the 
separations between the excited states. This behavior clearly departs from the Standard 2-level avoided 
crossing picture and is due to the contributions from the exponential number of terms in (^|) as will 
be analyzed elsewhere [|15| . We also note that the value n=10 does not correspond to the exponential 
scaling regime for the algorithmic complexity that appears to start for greater n vàlues as follows from 
the discussion above. 

In conclusion, we have performed numerical simulations of the adiabatic quantum optimization for 
SPP using a step-like density of states defined on a logarithmic scale of partition residues. The results 
indicate an exponential scaling of the algorithmic complexity as a function of the problem size. The 
apparent reason is the loss of structure in SPP during the effective coarse-graining over the intervals 
of partition residues corresponding to the same cost function vàlues. 
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